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It is well known that band insulators require an integer number of electrons (of each spin) per unit 
cell. Similarly, in bosonic insulators where the number of particles per unit cell (/) is fractional, the 
ground state either enlarges the unit cell, or realizes an exotic topologically ordered (fractionalized) 
state. Symmetric non-fractionalized Mott insulators only appear at integer /. However, the converse 
problem is relatively unexplored - is such a symmetric non-fractionalized insulator always allowed at 
integer /, or can it be prohibited by other lattice symmetries? An especially relevant example is the 
honeycomb lattice - where free spinless fermions at / = 1 (or, 1/2 site filling) are always metallic, 
due to point group symmetries. Here we argue that bosons at the same filling however can realize 
a Mott phase. We propose a wave function for this state and by a mapping to a classical partition 
function we compute its properties and demonstrate that the state is insulating, fully symmetric 
and has no topological order. Thus the absence of symmetry breaking in this case does not imply 
topological order. Our construction suggests that featureless insulators are generically allowed for 
bosons at unit filling on any symmorphic lattice in any dimension. 



Within the Landau paradigm, phases of matter are dis- 
tinguished by spontaneous symmetry breaking. Imphcit 
here is the assumption that a completely symmetric state 
(a paramagnet, in the context of spin systems) exists. 
In classical systems, this is just the high temperature 
phase. However, when discussing quantum phases at zero 
temperature, a paramagnetic state could be forbidden. 
For example, according to the Lieb-Schultz-MattiJi^ the- 
orem, for the spin-i Heisenberg antiferromagnet in one 
dimension, no completely symmetric quantum paramag- 
net exists. A spin disordered state with exponentially 
decaying spin correlations will necessarily break lattice 
symmetries. An extension of this theoreirP applies to 
two dimensions. In the square lattice spin-^ Heisenberg 
model, which has a half-odd-integer spin per unit cell, 
the spin disordered phase is not a trivial paramagnet. If 
it does not break lattice symmetries then it must be a 
quantum spin liquid, which has a hidden form of order 
called topological order. The latter leads to fractional- 
ized excitations with novel statistics, and is distinct from 
our notion of a featureless paramagnet. The absence of 
symmetry breaking can then be taken as indirect confir- 
mation of the quantum spin liquid, which is often used as 
a diagnostic both in numerics and in experiments- . Fur- 
thermore, in these systems where a trivial paramagnet 
is forbidden, quantum phase transitions often lie outside 
the Landau-Ginzburg- Wilson paradig nP. Therefore it is 
important to understand exactly when such trivial para- 
magnets are disallowed. 

These considerations can be readily translated to bo- 
son systems in a periodic lattic^, as realized by ultracold 
atomic gases in optical latticed. We assume a homoge- 
nous system, and the filling is defined as the number 
of bosons per unit cell. When the filling is not an in- 
teger, the ground state must break some symmetry, eg. 



by forming a superfluid or enlarging the unit cell, or, by 
realizing a topologically ordered state. On a simple lat- 
tice with say one site per unit cell, a Mott insulating 
state can appear at integer filling of bosonsP. This is 
the bosonic analog of the trivial paramagnet. In a sim- 
ple caricature of this state, each site is occupied by a 
fixed integer number of bosons. Clearly such a carica- 
ture does not work on more complicated lattices, with 
more than a single site per unit cell. For example, we 
can ask whether a lattice with n > 1 equivalent sites per 
unit cell can exhibit an insulating phase at a filling of 
one boson per unit cell. In some cases the lattice can be 
broken down into non-overlapping sets of finite number 
of sites, or molecules, which preserve all the symmetries, 
and can be considered as providing a single orbital at 
integer filling. However, for lattices like the kagome and 
honeycomb with n = 3, 2 respectively, no finite grouping 
of disconnected sites retains the symmetries of the lat- 
tice. In this case it is less obvious how to construct a 
Mott state at unit filling (which corresponds to 1/3 and 
1/2 boson per site respectively). One possibility is to 
consider a topologically ordered phase where emergent 
excitations carry a fraction of the boson charge, which 
can then be placed at integer site filling. However as we 
argue below, in these two situations and in several oth- 
ers, a fully symmetric non-fractionalized Mott insulator is 
possible. Throughout this paper, we consider only tight- 
binding models where particles are restricted to occupy 
sites on the given lattice. This is crucial to distinguish 
between the triangular, honeycomb and kagome lattices, 
which share the same space group symmetries. 

In a recent publication'' we demonstrated that a sym- 
metric and non-fractionalized Mott Insulator could be 
constructed on the kagome lattice at filling one per unit 
cell. This was accomplished by obtaining a family of ex- 
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ponentially localized Wannier orbitals that respect lattice 
symmetries, and occupying each one with a boson. This 
construction requires that a band insulator of spinless 
fermions can also be defined at the same filling, since a 
Slater determinant of these Wannier orbitals produces a 
band insulator. 

While a band insulator is possible at unit filling on the 
kagome lattice, it is known that (spinless) fermions on 
the honeycomb lattice are always gapless at filling one 
per unit cell (which corresponds to half a fermion per 
site). This is relevant to the band structure of graphene, 
which is gapless at the corners of the Brillouin zone (±K 
points), but also holds for more general tight binding 
Hamiltonians. The lattice symmetries of rotation and re- 
flection ensure that the two bands touch at these points. 
In this case our previous constructiorP is not applica- 
ble. In fact, one may suspect that the absence of a free 
fermion band insulator on this lattice may indicate that 
no fully symmetric bosonic insulator exists either. In this 
case, this would be a model which, even at integer unit 
cell filling, is prohibited from realizing a trivial insulating 
phase. 

However, we show that a separate wave function, which 
we term the Voronoi permanent, can be constructed on 
the honeycomb lattice and shown to have the desired 
properties i.e. it is an insulating wave function with the 
right filling, without either symmetry breaking or topo- 
logical order. Our key idea is simple - we construct a 
molecular orbital involving equal superposition of the six 
sites of the honeycomb, and insert one boson in each of 
these orbitals. The resulting wave function is a perma- 
nent of single particle wave functions that involve sites 
within the Voronoi unit cell (sites closest to the center 
of each hexagon), hence the term Voronoi permanent. 
Since sites are shared between hexagons, deducing prop- 
erties of this state requires further calculation. We show 
this wave function can be mapped to a two dimensional 
statistical mechanical problem of loops, whose proper- 
ties are calculated using Monte Carlo simulations. The 
result indicates that we have a fully symmetric and non- 
fractionalized Mott phase on the honeycomb lattice. 

Recently, a interesting numerical study of the honey- 
comb lattice Hubbard model with spin 1/2 fermions has 
appeared^. A phase at intermediate coupling appears 
to be insulating and spin gapped but without symmetry 
breaking of any kind, and was proposed to be a spin liq- 
uid with topological order. While this conclusion may 
ultimately be correct, and though the states we propose 
are not directly candidate wave functions with SU(2) spin 
symmetry, they are sufficiently similar that one may sus- 
pect that spin symmetric paramagnets may also exist on 
the honeycomb lattice, and may be an alternate expla- 
nation of these results. Positive signatures of a quantum 
spin liquid such as ground state degeneracy or topological 
entanglement entropjI^HS] should settle this question. 

We note that our method can be extended to any 
symmorphic lattice, where the point group symmetries 
and the translations are independently generated, to con- 



struct a candidate Mott state at unit filling. However 
for non- symmorphic lattices our method fails since we 
cannot identify a symmetric unit cell. At this point we 
may ask — is there any lattice on which a symmetric 
non-fractionalized Mott insulator is forbidden at inte- 
ger filling? While we do not have a definitive answer to 
this question, we note that non-symmorphic lattices are 
promising candidates. Familiar non-symmorphic lattices 
include the diamond and pyrochlore lattices. 

This paper is structured as follows. We begin by 
demonstrating the unique difficulties of constructing in- 
sulators on the honeycomb lattice, thus motivating our 
work. Next we explain our construction of the Voronoi 
permanent state, a candidate wavefunction for a feature- 
less bosonic insulator, and proceed to show that its cor- 
relations can be understood in terms of a classical loop 
model. We describe the simulation of this loop model us- 
ing a Monte Carlo worm algorithm, and use this to show 
that the state describes a symmetric insulating phase on 
the honeycomb lattice. We close with a discussion of the 
relevance of our results for quantum spin systems and 
extensions to other symmorphic lattices. 



I. FREE FERMIONS AND HONEYCOMB 
LATTICE SYMMETRIES 

Consider spinless fermions on the honeycomb lattice. 
The nearest-neighbor tight-binding model yields a spec- 
trum with Dirac points at the ±K points of the Brillouin 
zone, where two dispersing bands touch. At half-filling 
with one fermion per unit cell, the lower band is occu- 
pied and the Fermi energy is at the Dirac points. In 
order to gap out a Dirac point, a lattice symmetry (re- 
fiection, threefold rotation or inversion) must be spon- 
taneously broken. In fact the band touching at the ±K 
points holds beyond the Dirac limit, for any tight binding 
Hamiltonian on this lattice. To see this, we study the ir- 
reducible representations of the 'little group' of the Dirac 
points. At any momentum q in the Brillouin zone, the 
little group Gq is the subgroup of the space group that 
leaves q invariant or translates it by a reciprocal lattice 
vector; the Bloch Hamiltonian hq at q commutes with 
the little group generators and so we can classify energy 
bands using irreducible representations of Gq. 

Consider the K point; our arguments apply equally 
well to — K. Rotations by 2tt/3 leave K invariant, so 
the generator of threefold symmetry -R27r/3 = ^^/s 
the little group. The other little group generator is the 
mirror reflection a2 - which also leaves K invariant - and 
together these generate Gk = Ddh- These generators 
act on the two components of the Bloch function which 
represent the amplitudes to be on the two sublattices of 
the honeycomb. Their representation matrices take the 
form 
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These form a two-dimensional irreducible representation, 
and thus the band touching is protected by symmetry. 
This argument carries through for spinful fermions with 
SU(2) spin rotation symmetry in which case half-filling 
corresponds to two fcrmions per unit cell. Thus, there is 
no fermionic band insulator at half-filling on the honey- 
comb lattice that preserves all its symmetries. 

It was recently demonstrated that in some cases one 
can construct bosonic quantum Mott insulating wave- 
functions given the orthogonal, symmetric, and expo- 
nentially localized Wannier orbitals of a fermionic band 
insulator at the same filling. For the honeycomb, the 
protection of the band touchings by symmetry forbids 
the existence of Wannier orbitals that satisfy these crite- 
ria, and thus this approach fails. We therefore need an 
alternate construction which we discuss below. 



which preserves all point group symmetries. This is pos- 
sible on symmorphic lattices, for which all the symme- 
tries may be realized at a single point, as we elaborate be- 
low. For example, a similar hexagon state on the kagome 
lattice corresponds to 1/3 boson per site, since each unit 
cell has three sites. 

The expression ([T]) for |^'o) is more complicated than 
it would appear at first sight. Because hexagons on neigh- 
boring unit cells overlap, it is a highly entangled state. 
Such quantum entangled states have no classical ana- 
logue and arc difficult to visualize. More to the point, 
because of the nonorthogonality of the hexagons we can- 
not tell whether |^o) describes a featureless insulator, 
a superfluid or another symmetry broken state. This re- 
quires an explicit computation, to which we now turn. 



II. CANDIDATE HONEYCOMB MOTT STATE 

Consider the honeycomb lattice with one boson per 
unit cell. Since this corresponds to half a boson per site, 
we cannot construct a featureless insulator by restricting 
bosons to individual sites. The next option is to delocal- 
ize bosons across some supercell. However, any choice of 
non-overlapping supercells necessarily breaks symmetry. 
We could attempt to construct linear combinations of 
symmetric supercells that are orthogonal while preserv- 
ing the symmetry, but these would be symmetric Wan- 
nier orbitals which would not be exponentially localized 
because of the band touching on the honeycomb noted 
above. Instead, we use the smallest supercell that pre- 
serves the lattice symmetry, corresponding to the six sites 
of a hexagon. Neighboring hexagons have overlapping 
sites; this complication arises for any choice of supercell. 
This motivates the following 'hexagon product state', 
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Here j labels a site and R a unit cell, i.e. a site on 
the Bravais lattice. Or denotes the sites on hexagonal 
plaquette centered at R, and |0) is the boson vacuum. 
In first quantized form, the wavefunction is 



*o(i"i,r2,...,rAr) = perm[(/)R,(rj)] 



(3) 



where 



3r1L 
permanenD^ 



= (0|&r-SR|0) , and perm refers to the 
Note that since each hexagon is in one- 
to-one correspondence with a unit cell, and each unit cell 
has two sites, this state has the requisite 1/2 boson per 
site. Since they are permanents of bosonic orbital^^^ de- 
fined on a symmetric Voronoi cell of the lattice, we term 
such wavefunctions Voronoi permanents. 

The key to this Voronoi permanent construction is 
to associate with each unit cell a single-particle orbital 



III. MAPPING |*o) TO A LOOP MODEL 

We wish to study properties of |^'o) such as the boson 
correlation function: 



G''°^(z,j) = 



(vl/o|&l&,|*o) 



(4) 
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which describes the ability of the particles to propagate 
between lattice sites i and j. In a Mott insulator, this 
decays exponentially, while for a 2D superfluid at T = 
it exhibits long range order. Actually as we demonstrate 
below we can map the correlations in |V-'o) to those of a 
classical finite temperature 2D model with short-range 
interactions. True long-range order is thus ruled out 
by the Mermin- Wagner theorem. Correlations must de- 
cay at least algebraically, corresponding to a Kosterlitz- 
Thouless (KT) superfiuid phase. The only other possible 
ordering is a discrete breaking of lattice symmetry, which 
can be diagnosed by studying the spatial structure of cor- 
relations. Topological order is ruled out for this state 
since it is an unconstrained product state. 

Relating the ground-state wavefunction to Boltzmann 
weights of a classical statistical mechanical model has no- 
table precedents including the Laughlin fractional quan- 
tum Hall wavefunctioiP^, the Rokhsar-Kivelson wave 
function of dimer modelJ^^E^ and the AKLT spin 
wavefunctiorP^ll^ni^ In our case there are different ways to 
represent G*"* in terms of a classical partition function. 
The first involves the statistical mechanics of closed loops 
on the triangular lattice, and since it offers multiple ad- 
vantages it will be our primary tool. A second approach 
(described in the supplementary material) is based on co- 
herent states of bosons, and was used to benchmark loop 
model results. 

We begin by showing that the normalization (^'ol^o) 
is the partition function of the loop model derived below. 
Following the definition of |^o) we have 

(^Po|*o) = (0| (Dr^r) I (Dr' |0), (5) 



where the products are over sites R of the (triangular) 
Bravais lattice. By commuting Br, Bjj^ past each other. 
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we find that the correlation function is mapped to a sum 
over products of commutators in which each i?R operator 
is matched with one -BJj^ operator, 



(*oi*o>=En[^i^'<R)] 



(6) 



R 



where a denotes a permutation of the sites on the tri- 
angular lattice. Since only neighboring hexagons share 
sites, the commutator is 



-Br, S^,] 5r,R' + TO(5R/^„n(R) 



(7) 



where nn(R) denotes the nearest neighbors of R on the 
Bravais lattice. Here m,— \ (m = i) for the honeycomb 
(kagome) lattice. In general, m = "pjq where p is the 
number of sites shared by a pair of neighboring Voronoi 
cells, each of which has q sites. For a hexagonal Voronoi 
cell g = 6, and p ~ 2(1) for the honeycomb (kagome) 
lattice. 

Using I?]) we can restrict the sum over permutations 
in (|6]) to those in which each site is matched with either 
(i) itself, contributing a multiplicative factor of 1 to the 
weight (which we represent as an empty site), or (ii) a 
neighboring site, contributing a factor m (which we rep- 
resent as an arrow pointing from B to -B^). Since every 
site must be matched to exactly one other site, the ar- 
rows form closed loops which cannot intersect or touch; 
a loop of length L contributes a multiplicative factor of 
to the partition sum. Thus, the normalization defines 
a statistical mechanical model of closed, nonintersecting 
directed loops (which we take to include empty sites and 
dimers) on the triangular lattice, in which each loop seg- 
ment contributes a multiplicative weight m to the prob- 
ability of a given configuration. The inclusion of empty 
sites and dimers, and the constraint that loops cannot in- 
tersect (or touch), distinguish this particular loop model 
from more conventional ones studied. 

A similar argument relates the numerator of Q to a 
defect correlator of this loop model: the and b opera- 
tors each remove one of the B and B^ operators respec- 
tively, giving a pair of defects, a site with an arrow out 
but no arrow in and its converse. Thus each configura- 
tion contributing to the numerator must include a single 
open loop. 

This interpretation as a loop model unifies the states 
on all lattices with the same underlying Bravais lattice. 
For the triangular Bravais lattice, the physical states on 
the kagome and honeycomb appear as specific points in 
the continuum of possible values for the loop weight m. 
For an intuitive understanding of the meaning of this 
continuum of values, consider decorating the honeycomb 
lattice by adding a site (represented by a^) to the center 
of each hexagon, and modify the definition of bI^, 
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R 



cos 9a 



sine . 
V6 ^ 



(8) 



Since only b sites are shared by hexagons, computing 
the commutator now yields ilj with m = | sin^ 6 on the 




FIG. 1; Voronoi construction and loop mapping 

(L) Decorated honeycomb and delocalization of bosons onto 
hexagons as in Eq. ([8|; m increases from to 1/3 from top 

to bottom. (R) Sample loop configuration showing loops, 

dimers and empty sites which do not touch or intersect. 

Below we show the continuum of loop weights m, where 
m = is the atomic insulator and m = 1/6, 1/3 correspond 

to Voronoi permanents on the kagome and honeycomb. 



honeycomb. For = 0, we have an atomic insulator, 
since the bosons are restricted to the central site. As 
we increase 9 we spread bosons across hexagons until we 
arrive at the 'honeycomb point' 9 — 7r/2 (to — 1/3) where 
we can remove the empty central site. We have sketched 
this construction of the Voronoi permanent in Fig. [l] 

Thus, the mapping to the loop model lets us inter- 
polate from the atomic (m = 0) limit to the honey- 
comb (to — 1/3) point as shown schematically in Fig. T] 
This interpolation allows us to study the evolution of the 
correlation functions from the atomic insulator limit at 
TO = where the correlation length vanishes. 

We compute correlation functions of the loop model us- 
ing a modified classical Monte Carlo worm algorithnPH. 
The worm algorithm is especially suited to this problem 
because it works naturally with the loop representation 
and as noted previously, breaking continuous symmetry 
in the model yields a KT superfluid with algebraically 
decaying correlations. First, unlike most approaches the 
worm algorithm is not expected to exhibit critical slow- 
ing down in the algebraic phase. Second, it can often 
access large system sizes more efficiently than alterna- 
tive approaches. Third, knowledge of the worm winding 
numbers allows us to directly compute the helicity mod- 
ulus (proportional to the superfluid density) , a definitive 
diagnostic of the KT phase^. This vanishes in the disor- 
dered phase and exhibits a universal jump of ^ at the KT 
transitiorP^l. More complicated superfluid orders (e.g. a 
'pair superfluid' where ((6-)^) ^ but (6|) = 0) wiU also 
be captured in this approach. 

We performed extensive worm algorithm Monte Carlo 
simulations for periodic L x L triangular lattice systems, 
with L = 12, 16, 20, and averaged over 10^ Monte Carlo 
steps per site (MCS) in each case (10^ MCS were suffi- 
cient for the kagome at low to = 1/6). Each system was 
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FIG. 2: Loop model correlation function. 

Linear-log plot with distance measured along a basis vector. 
The fastest possible KT algebraic decay, ~ r~^/*, is shown 
for comparison; m = 1/3 (corresponding to the honeycomb 
lattice) displays exponential decay, indicating an insulating 
phase. Error bars are smaller than the line widths. 

initialized in the 'm = 0' configuration with no loops, and 
allowed to equilibrate over 50 000 MCS before record- 
ing averages. Error estimates were obtained by comput- 
ing the standard error of the data from 20 independent 
runs, a conservative approach. We benchmarked results 
against coherent state simulations at m = 1/3,1/6 and 
loop perturbation theory at small rrPH and found excel- 
lent agreement. Our results are summarized in Figs. [2j 
landH 



IV. NUMERICAL RESULTS 

We study the properties of the honeycomb state in two 
distinct ways. First, we focus directly on the honeycomb 
lattice. In addition to computing the correlation func- 
tion and winding numbers (hence, the superfluid density) 
in the loop representation using the worm algorithm (at 
m = 1/3), we also compute the correlation function for 
smaller system sizes directly using a coherent state repre- 
sentation of the honeycomb wavefunction (see the supple- 
ment for details). In addition, we determine whether any 
discrete symmetry is broken by explicitly checking that 
the correlations preserve all lattice symmetrieJ^Sl Second, 
we enlarge the study to include other points on the m 
axis, including to = 1/3 for the honeycomb, to = 1/6 
for the kagome, to = 1/4 as an intermediate state and 
TO = 17/48 which lies beyond the honeycomb point. We 
use the worm algorithm to study the evolution of the 
superfluid density and the correlation function from the 
trivial to = atomic insulator to the to = 1/3 honey- 
comb and even slightly beyond. 

We find that G^°''{i,j) decays exponentially with a cor- 
relation length of ^ ~ 2.4 lattice sites, as shown in Fig. 
[2j The correlation length evolves from ^ = at to = 
to ^ ^ 0.9 lattice sites at the kagome (to = 1/6) and 
remains a small fraction of the system size for to values 
beyond the honeycomb. As is clear from Fig. [2] the 



FIG. 3: Helicity modulus. 

The helicity modulus (proportional to the superfluid 
density) tends to zero for increasing L, consistent with the 
expectation from KT finite-size scaling for the insulating 
phase. (Inset) Same figure, with axes zoomed out to show 
that the helicity moduli for all m values studied are much 
smaller than the universal value of 2/-K at the KT transition 
(the smallest superfiuid density allowed in the KT phase). 
Error bars are smaller than the markers. 

correlations decay exponentially and cannot be fit by a 
power law. We also plot for comparison the most rapid 
possible algebraic decay in the KT phase, G(r) ^ r^^/"*, 
to emphasize that algebraic decay of single boson corre- 
lations is ruled out from this analysis. 

As shown in Fig. |3j the superfluid densities are all 
much less (see inset) than the universal jump value 2/7r 
at the KT transition^ indicating that the wavefunction 
remains in the insulating phase for the to values studied. 
Multi-boson condensates are also ruled out by this ob- 
servation. As expected, we find that smaller system sizes 
overestimate the superfluid order, and the superfluid den- 
sity decreases with increasing size, consistent with the 
claim that the honeycomb Voronoi permanent describes 
an insulating state. 

Finally, neither the loop model nor the coherent state 
simulations exhibit breaking of discrete lattice symme- 
tries for the honeycomb lattice Voronoi permanent. To 
determine this, short MC runs were used to avoid aver- 
aging out any symmetry-breaking by oversampling. As 
an example, a 500 MCS coherent boson run for L = 12 
gave rms values of symmetry-breaking of less than 2% 
of the maximum. As a simple visual demonstration of 
lattice symmetry, a sample correlation function for the 
loop model on the triangular lattice at to = 1/3 from an 
L — 12, 500 MCS worm algorithm run is shown in Fig.|4] 
Details of the analysis are provided in the supplementary 
material. 



V. EXTENSIONS TO OTHER LATTICES 

Our construction of a candidate symmetric Mott insu- 
lator wave function at unit filling can be extended to any 
symmorphic lattice, where a primitive unit cell that re- 
spects all lattice symmetries is present (the Voronoi cell 



FIG. 4: Lattice symmetries in the correlation 
function. 

Interpolated contour plot of loop model correlation function 
at m = 1/3, with white space in the center corresponding to 
the central peak. Correlations decay rapidly and are 
consistent with lattice symmetries. 



for example). Imagine introducing a site at the center 
of the unit cell and filling it with one boson (the atomic 
insulator). Then, gradually transfer the boson to the 
physical sites in the unit cell thereby obtaining a family of 
wave functions controlled by one parameter as in Eq. ([8| . 
This family of wavefunctions has a loop model descrip- 
tion defined on the underlying Bravais lattice. When the 
added sites are fully depleted they can be detached to 
yield the desired Voronoi permanent wave function. The 
loop model description of the state can then be used to 
check whether the Voronoi permanent lies in the same 
phase as the atomic insulator, using techniques similar 
to those in this paper. Symmorphic lattices in any di- 
mension thus always admit the construction of Voronoi 
permanents. Determining if the state is ordered must be 
done on a case-by-case basis. Here we described the most 
basic construction — more complicated wave functions 
with additional tuning parameters could be constructed 
on similar lines, which could expand the range of avail- 
able states. 

Another 2D example is the checkerboard lattice with 
two sites in the unit cell, which may be viewed as a set 
of corner sharing tetrahedra in two dimensions. Again, 
the (spinless) fermonic bands must touch due to lattice 
symmetriea^Sl. However, our Voronoi permanent con- 
struction could be applied here too. An example in d = 3 
is the cubic perovskite lattice of corner-sharing octahedra 
( Fig[5| with three sites per unit cell. Here fermions are 
again gapless at filling one per unit cell but the Voronoi 
permanent construction applies here as well. Another 2D 
example with a band touching was considered in Ref. 1271 
but in this case it is trivial to construct a Bose insulator 
by filling non-overlapping sets of sites. 

On non-symmorphic lattices it is impossible to choose 
a unit cell that respects all the space group symmetries 
(or alternatively, to choose a unique symmetric point at 
which to place the decorated site) and the construction 
fails. Examples of non-symmorphic lattices include the 
well known pyrochlore and diamond lattices. 



FIG. 5: Corner sharing octahedron lattice 

This is a symmorphic lattice with sites at the vertices of 

corner-sharing octahedra. A Voronoi permanent 
wavefunction of a candidate insulating state at unit cell 
filling of one (corresponding to 1/3 site filling) can be 
constructed. 

VI. DISCUSSION 

We began this paper by posing a question: is a fea- 
tureless and non-fractionalized (i.e., topologically triv- 
ial) insulating phase possible for a system of bosons on 
the honeycomb lattice at a filling of one boson per unit 
cell? By explicitly constructing a simple trial wavefunc- 
tion and computing correlations within it, we have ar- 
gued that this is indeed the case: the Voronoi permanent 
wavefunction we proposed captures the essential prop- 
erties of such a featureless insulating phase. The hon- 
eycomb is an especially interesting example, because its 
free fermion band structure cannot be insulating with- 
out breaking symmetry. As a corollary, we found that a 
similar wavefunction on the kagome lattice — which can 
be viewed as the truncation of wavefunctions of Ref. [7] 
— also represents an insulating phase. If we attempt 
to make a Slater determinant of the honeycomb lattice 
Voronoi orbitals, the state vanishes. However, a similar 
procedure on the kagome lattice is expected to lead to 
a band insulating state of fermions which is allowed by 
symmetry. 

In both cases the wave functions obtained are positive 
definite, and could be ground states of an unfrustrated 
bosonic model. Here, we have not discussed Hamilto- 
nians for which these are ground states. However, gen- 
eral methods to construct parent Hamiltonians for wave 
functions like those discussed here — which are cap- 
tured by partition functions of local statistical mechan- 
ical models — have been discussed by HenlejP^. Fur- 
thermore, recently ultracold atoms have been confined 
to kagom^^H and honeycomb?- optical lattices. This mo- 
tivates the study of realistic Mott-Hubbard model Hamil- 
tonians which could potentially realize such fractional 
site filling insulators. Both these directions are left to 
future work. 

Our wave function can also be interpreted as a spin 
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wave function where different boson occupations are dif- 
ferent Sz values, and only U(l) spin rotation about the 
z axis is a symmetry. The honeycomb lattice wave func- 
tion then represents a featureless magnetization plateau 
at 2/3 of the saturation value, with S — 3/2 spins on 
the sites. We can also write down a wave function for 
S = 1/2 by projecting out multiple boson occupation 
from our state (hardcore limit). This state is appropri- 
ate for a U{1) symmetric spin system, and, if projec- 
tion does not modify its basic properties, it represents 
a spin gapped Sz — state on the honeycomb lattice 
without topological order. The effect of projection must 
be checked numerically, but is expected to further sup- 
press superfluid correlations and not affect the topologi- 
cal properties, although it may enhance lattice symmetry 
breaking effects. 

Our approach is readily extended to all symmorphic 
lattices, where analogues of the hexagon product state 
can be immediately constructed. While numerical checks 
are required to confirm these are states with the desired 
properties, these are straightforward since they reduce 
to classical statistical mechanics of loops. We note that 
it is not even possible to construct such promising can- 
didates for bosons at fractional unit cell filling. This 
suggests that the unit filling of symmorphic lattices is 
an essentially different situation, where symmetric non- 
fractionalized Mott states are generically allowed. 

We note however that out method fails for non- 
symmorphic lattices, which include interesting examples 
such as the diamond and pyrochlore lattices. We spec- 
ulate that a fully symmetric non-fractionalized bosonic 
Mott state is absent on such lattices, even at a filling of 
one boson per unit cell. 
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SUPPLEMENTARY MATERIAL 

1. Loop Model Mapping 

Let us recapitulate the mapping from the quantum bo- 
son correlator to the statistical mechanical loop model in 
more detail. For both the kagome and honeycomb lat- 
tices, it is convenient to explicitly decompose the coordi- 
nate i = (R, a) where a is a sublattice index and R ~ 
msLi + na.2 is a point on the triangular Bravais lattice. 



with ai and a2 primitive lattice vectors (shown in the fig- 
ure) and we define for future reference = — (ai -f a2). 
We wish to compute the boson correlation function 

G^-,(R,R') = (&La&R',c') (9) 
within the hexagon product state/Voronoi permanent 

R \R',a J 

For the kagome lattice, a = 1,2,3 where we choose a 
symmetric unit cell in which the three sites belonging to 
a unit cell lie in the directions of the primitive Bravais 
lattice vectors ai,a2,a3 respectively. For the hexagon 
state we then have /R(R',a) = ^ ((5r3' -l-(5R3'+e„) 
where the final refers to the Bravais lattice vector 
pointing from one hexagon to an adjacent hexagon which 
shares sublattice site a with it. For our choice of unit cell, 
— a^ . 

For the honeycomb, a = 1,2 where the unit cell 
is taken to consist of the the two sites on a vertical 
bond. For the hexagon state we then have /r(R',q;) — 
^ (<5r,r' + fe^R'-ai + ^R,R'+v„) where we define vi = 
a3 and V2 — &2- 

Again we begin with the normalization ( i.e. the de- 
nominator of Q) 

(*o|*o) - (0| (UB^'j (UbI^ |0) (11) 

withi?t = ER',a/R(R':«)6k',a- 

A single pair B^^^, , i?R gives the commutator 
C[R,R'] = ER,3..oi,a./i^(I^i'"i)/R'(I^2,a2). In a 
Wick type decomposition, this becomes a sum over all bi- 
jective (one-to-one) maps i.e. permutations a : R — )■ R', 
with each term in the sum being the product of com- 
mutators Yiu, ^[£1, o'(R)]- We are saved from computing 
this functional integral because the hexagons of two unit 
cells overlap only if they belong to the same or neighbor- 
ing unit cells, C[R, R'] = (5r,r' + TO(5r/ „„[r], so maps 
a only appear in the sum if they take R either to itself 
(weight 1) or to a neighbor (weight m). 

Making contact with the description in the main text, 
an allowed map a can be pictured as a collection of ar- 
rows between neighboring sites on the triangle lattice, 
with each site R either having no arrows (in case a maps 
it to itself) or exactly one arrow going out (to o'(R)) and 
one arrow coming in (to tT~^(R)). Therefore (^'ol^'o) is 
a sum over directed closed loop configurations on the tri- 
angle lattice, with each configuration in the sum weighted 
by m to the power of the total length of its loops. Loops 
may not touch or intersect with the exception of the zero 
area length two loop, which is permitted. Strictly speak- 
ing this is gives generalized vertex model (with 37 states 
per vertex, 6^ two-in-two-out states and one for the case 
when the site is mapped to itself) and not a simple loop 
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model, though we may stih study it by the worm algo- 
rithm by incorporating two bond-occupation flavors. We 
are interested in m = 1/6 and m = 1/3, corresponding 
to the kagome and honeycomb lattices respectively. 

Next, we turn to the numerator. Each b operator 
knocks out a B operator with a coefhcient /, resulting in 
a sum over defect correlators, 

Gi%,{R,R')= J2 /R,(R,«)/R.(R',a')G™(Ri,R2) 

R,i ,R,2 

(12) 

where we named the defect correlator G^°™ because of 
its natural interpretation as a correlation function in the 
worm algorithm for Ri ^ R2, 

(13) 

For the kagome, the boson correlator is given by the sum 
of four worm correlator values; for the honeycomb, the 
sum of nine. Thus, we can extract correlation functions 
on the honeycomb and kagome lattices explicitly from 
the worm algorithm defect correlator. 

2. Loop Perturbation Theory 

Within first-order perturbation theory in m, loops are 
costly and the sum in the numerator will involve only 
configurations where the defects are connected by the 
shortest line segment between them. To this order, re- 
stricting for a moment to a ID lattice, we find that de- 
fect configurations are all multiplied by an extra factor of 
^\R-R I ^ xhis naive estimate holds when there is only one 
shortest line connecting the defects, as is the case when 
R — R' is a multiple of a triangle lattice vector i?i or i?2. 
However, consider the case where R — R! ^ n{Ri + i?2)- 
Then there are many shortest line segments between the 
defects, identical to the number of shortest paths between 
opposite corners of an n by n square lattice. There are 
(^^) such paths, asymptotically going as 

4|i?,-fl'|_ Thus 

we find that the asymptotic behavior of the correlation 
function is 

(blh) ^ {Amf 

within first order perturbation theory in m ^ l/T. This 
perturbation theory result turns out to overestimate the 
extent of the algebraically decaying superfiuid phase. 
However, it provided a further benchmark of the worm al- 
gorithm computation of the loop model correlation func- 
tion. 



3. Worm algorithm calculations 

The directed loop configurations with and without an 
open loop may be studied simultaneously using a Monte 



Carlo 'worm' algorithrrPU . The basic idea of the worm al- 
gorithm is to simultaneously gather statistics on the cor- 
relation function and the normalization in Q by work- 
ing directly in the loop representation. The 'worm' is an 
open loop configuration; by allowing a worm to shrink or 
grow by a random process, until it closes of its own ac- 
cord while preserving detailed balance - implemented by 
a local-update Metropolis algorithm - one gains statis- 
tics on both open and closed configurations, which con- 
tribute respectively to the numerator and denominator 
of Q. The canonical example of the use of a worm al- 
gorithm is to study the XY model in the loop-current 
(dual) representation. 

Unlike the usual XY case, our loop model has strong 
interactions in that the loops are forbidden to touch or 
intersect, necessitating some modifications to the algo- 
rithm. Here, we merely note that the twin complications 
of loop self-avoidance and the triangular lattice geom- 
etry presents unique challenges for the worm algorithm. 
Specifically, at higher values of m when long loops are rel- 
atively favorable, the worm can occasionally get 'stuck' 
in a configuration (for instance, a spiral) for which nearly 
all proposed updates will be rejected. To extricate the 
worm from such a configuration requires exponentially 
many update steps and thus the correlation function is 
poorly sampled and the algorithm fails to converge. This 
issue necessitated long runs ( > 10® MCS) for the hon- 
eycomb, though for smaller m we found 10^ MCS to be 
sufficient. For the honeycomb (m — 1/3) a single run at 
the largest system size (L = 20) did not converge even af- 
ter 10^ updates; we discarded it and a similar m = 17/48 
run when taking averages. We note that the honeycomb 
data for L = 12, 16 are free of such convergence issues. 

4. Additional Numerical Data 

In Fig. |6] we show a log-log plot of the worm algorithm 
correlation function to a power law for the largest system 
size, L = 20 (same data as in Fig. [2| with the fastest KT 
algebraic decay shown for comparison on a log- log plot. 
It is clear that a power-law fit is inconsistent with the 
data as this would appear as a straight line on this plot. 

5. Helicity Modulus in Terms of Winding Numbers 

The relation between worm winding numbers and the 
helicity modulus may be seen as follows. Take the loop 
model on a torus and thread a fiux 9 through one of 
its handles, say the cycle associated with the periodic 
boundary condition r ~ r -I- Lx in the x direction. This 
is equivalent to a uniform vector potential Ax perme- 
ating the system with magnitude A = 0/L. The Boltz- 
mann factor for the worm to grow a step 5R now appears 
multiplied by the phase exp(zy4.i; • 5R). This phase fac- 
tor cancels out for any closed loop, unless it threads the 
torus, crossing the periodic boundary conditions with a 
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states that &i|{zi}) = Zi|{zi}) and (0|{zi}) = 1 we find 



SF, T = Tkt 



honeycomb 



kagome 



1.5 2.0 



5.0 7.0 



FIG. 6: Loop model correlation function. 

Log-log plot as a function of distance along a basis vector, 
demonstrating power-law fit is impossible. Error bars 
smaller than linewidths. 



nonzero winding number Wx = J x ■ SR/L. The winding 
number is integer-valued for closed loops. Let Wx stand 
for the sum of all x winding numbers in a given closed 
loop configuration. Then the free energy owing to the 
flux, AF[e] = F[e]- F[0] is given by 



(bib,) = Sjo + 



J [dz,dz-,]e-^^\^'\' UTL\l3n{z)\Hz,-zo) 



J [dzidz 



'E.N.I" 



(18) 



where Bu.{z) — X^i /R.(*)-^i- We see that it suffices to 
consider the correlation function (zrZq) of the statisti- 
cal mechanical model described by a complex number at 
each site, with the the manifestly positive semi-definite 
Boltzmann weight r{{z,}) = Hr \Bn.iz)\'^e-^i l^'l'. 

As the model has only U{1) symmetry (as is evident 
from the form of the wavefunction) , even in two dimen- 
sions it can exhibit an algebraically correlated KT phase, 
and so the worm algorithm Monte Carlo on the loop 
model is a more convenient tool to study the phase. Nev- 
ertheless, the coherent boson approach can be readily 
simulated using a simple local-update Metropolis algo- 
rithm which allows us to benchmark the worm results on 
small system sizes. We find excellent agreement between 
the worm and coherent boson results through the largest 
system size studied with the latter algorithm, L = 12. 



exp{-/3AF[e]) = {exp{i9Wx) 



(14) 



7. Discrete Symmetry Breaking 



yielding the helicity modulus 

Y = d^F/d9^\e^o = {{W) - (W)^) . (15) 

The superfiuid inverse temperature /? which appears in 
these expressions is set to 1 in our model. Studying wind- 
ing number fluctuations (or, the combination /3Y) thus 
allows us to identify the transition point. 

6. Bosonic Coherent States 

Another way to map the wavefunction |vI/o) is to work 
in the basis of bosonic coherent states. Let us consider 
the bosonic correlation function, 

G-(o,i).(.S.,).., + ff2M^ ,10) 

We now use the resolution of identity for bosonic coherent 
states, 



[dz,dz-,]e-^»l-'l'|{za)({z,}| (17) 



dl{.czidlmzi 



In- 



where = and [dz^dzi] = JJ^ 

serting this above, and using the property of the coherent 



The worm algorithm explicitly computes the correla- 
tion functions of the boson order parameter, and we are 
able to conclude that there is no long-range superfluid 
order. However, this does not immediately rule out dis- 
crete symmetry-breaking which could occur via an Ising 
transition. To study this from the correlation function, 
it is convenient to work with the correlation function 
Gl°^{R) = {bli ,^boj3) from the site 0,/3 to the site R, a 
where as before we label the sites by their Bravais lattice 
vector and sublattice index. If we now perform a discrete 
(i.e. point-group) symmetry operation U on the lattice, 
we have in general that G^°|(R) ^ UG''^^(R)U-^ with 
the operator U transforming both the site index as well 
as possibly the sublattice indices. We may deflne the 
'Ising' order parameter corresponding to symmetry U as 



(r; 



(19) 



and by explicitly examining the correlation functions we 
can verify that this vanishes asymptotically within our 
error estimates. In the main text we quote the rms aver- 
age of X^^-* (R). We perform this test explicitly for mirror 
reflections and rotations that together generate the lat- 
tice point group. For these symmetry breaking tests we 
use short Markov chains of 100 - 2000 MCS, to ensure 
that the chains do not sample multiple symmetry broken 
states while still collecting sufficient statistics for comput- 
ing the order parameter. To check for symmetry break- 
ing we primarily used simulations of the wavefunctions 
in the coherent state representation, since fluctuations in 
the associated Markov chains exhibit shorter correlations 
than those in the worm algorithm. 
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